home *** CD-ROM | disk | FTP | other *** search
/ Shareware Grab Bag / Shareware Grab Bag.iso / 090 / contour.arc / CONTOUR.BAS
Encoding:
BASIC Source File  |  1987-06-11  |  7.2 KB  |  199 lines
  1. 1 '*** Contour Map example: CONTOUR1.BAS
  2. 2  CLS: PRINT "CONREC Example. Coutour the function ":PRINT
  3. 3  PRINT "f(x,y)=sin((x^2+y^2)^0.5) + ((x-c)^2+y^2)^-0.5 ":PRINT
  4. 4  PRINT "letting x and y range from -2pi to +2pi.":PRINT
  5. 5  PRINT "Building the data set now.  Please wait..."
  6. 6  OPTION BASE 0                 'lower bound of zero for all array indices, instead of 1
  7. 7  PI#=3.141592654#
  8. 8  TRUE=-1: FALSE=0
  9. 9  ILENGTH=319: JLENGTH=199      'dimensions of output plot axes (CGA full screen)
  10. 10 IMIN=0: JMIN=199             'coords of left bottom corner of screen
  11. 11 IUB=30: JUB=30: NC=10        'number of grid intervals and contour levels
  12. 12 DIM D(IUB,JUB),X(IUB),Y(JUB) 'data array
  13. 13 DIM Z(NC-1)
  14. 14 ' define function and coords
  15. 15 FOR I=0 TO IUB               'check all x-grid levels
  16. 16  IX=2*PI#*(2*I-IUB)/IUB       'ix ranges from -2pi to +2pi
  17. 17  FOR J=0 TO JUB              'check all y-grid levels
  18. 18    JY=2*PI#*(2*J-JUB)/JUB     'jy ranges from -2pi to +2pi
  19. 19    R=SQR(IX^2+JY^2)
  20. 20    D(I,J)=SIN(R)+0.5/SQR((IX+3.05)^2 + IY^2) 'evaluate user supplied function,fill array
  21. 21  NEXT J
  22. 22  X(I)=I*ILENGTH/IUB+IMIN     'scale x(i) to span plot area
  23. 23 NEXT I
  24. 24 FOR J=0 TO JUB
  25. 25   Y(J)=JMIN-J*JLENGTH/JUB    'scale y(i) to span plot area
  26. 26 NEXT J
  27. 27 FOR I=0 TO NC-1
  28. 28   Z(I)=(I-5)/5               'contour levels at -1,-0.8,-0.6, ... ,1
  29. 29 NEXT I
  30. 30 '
  31. 31 CLS: SCREEN 1,0              'set CGA screen 320 x 200
  32. 32 '
  33. 33 LINE(IMIN,JMIN-JLENGTH)-(IMIN+ILENGTH,JMIN),,B    'use a box for axes
  34. 38 GOSUB 100                    'Do conrec subroutine
  35. 40 IF NOT(PRMERR) THEN PRINT:PRINT:PRINT MSG$;
  36. 45 PRINT CHR$(7);               'alert user that program has finished output
  37. 50 WHILE LEN(INKEY$)=0: WEND    'press any key to finish
  38. 60 END
  39. 70 '
  40. 80 '
  41. 99 '
  42. 100 '**************************** CONREC.BAS *******************************
  43. 101 'Source:  Translated from original program by Paul Bourke,
  44. 102 '         "A Contouring Subroutine", p.145, BYTE Magazine (June 1987)
  45. 103 'Input:
  46. 104 '         d(0:iub, 0:jub)     matrix for data surface
  47. 105 '         iub,jub             index bounds of data array
  48. 106 '         x(0:iub)            array for column coords
  49. 108 '         y(0:jub)            array for row coords
  50. 109 '         nc                  number of contour levels
  51. 110 '         z(0:nc-1)           contour levels to use, in increasing order
  52. 111 '         VECOUT              external subroutine to plot contour lines
  53. 112 '
  54. 113 '
  55. 180 FALSE=0: TRUE=-1
  56. 200 DIM H(4)                  'releative heights of box above contour
  57. 201 DIM ISH(4)                'sign of h()
  58. 202 DIM XH(4)                 'x coords of box
  59. 203 DIM YH(4)                 'y coords of box
  60. 204 DIM IM(3)                 'mapping from vertex numbers to x offsets
  61. 205 '
  62. 206 IM(0)=0: IM(1)=1: IM(2)=1: IM(3)=0
  63. 207 DIM jm(3)                 'mapping from vertex numbers to y offsets
  64. 208 JM(0)=0: JM(1)=0: JM(2)=1: JM(3)=1
  65. 209 DIM CASTAB(2,2,2)         'case switch table
  66. 210 '
  67. 220 DATA 0,0,8,0,2,5,7,6,9,0,3,4,1,3,1,4,3,0,9,6,7,5,2,0,8,0,0
  68. 221 '
  69. 230 FOR K=O TO 2: FOR J=0 TO 2: FOR I=0 TO 2
  70. 240   READ CASTAB(K,J,I)
  71. 250 NEXT I: NEXT J: NEXT K
  72. 251 '
  73. 253 'Validate input params (PRMERR is parameter error flag)
  74. 254 '
  75. 260 PRMERR=FALSE
  76. 270 IF (IUB<=0)OR(JUB<=0) THEN PRMERR=TRUE
  77. 280 IF NC<=0 THEN PMERR=TRUE
  78. 290 FOR K=1 TO NC-1
  79. 292   IF (Z(K)<=Z(K-1)) THEN PRMERR=TRUE
  80. 294 NEXT K
  81. 295 IF PRMERR THEN MSG$="Error in input parameters":RETURN
  82. 297 '
  83. 298 'Scan array, top down and left to right
  84. 299 '
  85. 300 FOR J=JUB-1 TO 0 STEP -1     'downto
  86. 310  FOR I=0 TO IUB-1
  87. 311  ' find lowest vertex of the four
  88. 320  IF (D(I,J)<D(I,J+1)) THEN DMIN=D(I,J) ELSE DMIN=D(I,J+1)
  89. 321  IF (D(I+1,J)<DMIN) THEN DMIN=D(I+1,J)
  90. 322  IF (D(I+1,J+1)<DMIN) THEN DMIN=D(I+1,J+1)
  91. 324  ' find the highest vertex
  92. 325  IF (D(I,J)>D(I,J+1)) THEN DMAX=D(I,J) ELSE DMAX=D(I,J+1)
  93. 326  IF (D(I+1,J)>DMAX) THEN DMAX=D(I+1,J)
  94. 327  IF (D(I+1,J+1)>DMAX) THEN DMAX=D(I+1,J+1)
  95. 328  '
  96. 329  IF (DMAX<Z(0) OR DMIN>Z(NC-1)) THEN GOTO 9000
  97. 330  '
  98. 331  'Draw each contour withing this box
  99. 332  FOR K=0 TO NC-1
  100. 340    IF (Z(K)<DMIN) OR (Z(K)>DMAX) THEN GOTO 8000
  101. 350    FOR M=4 TO 0 STEP -1
  102. 360      IF M>0 THEN GOSUB 11000
  103. 370      IF M=0 THEN GOSUB 12000
  104. 380      IF H(M)>0 THEN ISH(M)=2 :ELSE IF (H(M)<0) THEN ISH(M)=0 :ELSE ISH(M)=1
  105. 390    NEXT M
  106. 391 '
  107. 400 '  Scan each triangle within box
  108. 401 '
  109. 410    FOR M=1 TO 4
  110. 420      M1=M: M2=0: M3=M+1: IF (M3=5) THEN M3=1
  111. 425 '
  112. 430      CASE = INT(CASTAB( ISH(M1),ISH(M2),ISH(M3) ))
  113. 440      IF CASE=0 THEN GOTO 7000
  114. 450      ON CASE GOTO 1100,1200,1300,1400,1500,1600,1700,1800,1900
  115. 460 '
  116. 1100     'line between vertices m1 and m2 (case 1)
  117. 1110     X1=XH(M1): Y1=YH(M1): X2=XH(M2): Y2=YH(M2)
  118. 1120     GOTO 6000
  119. 1130 '
  120. 1200     'line between vertices m2 and m3 (case 2)
  121. 1210     X1=XH(M2): Y1=YH(M2): X2=XH(M3): Y2=YH(M3)
  122. 1220     GOTO 6000
  123. 1230 '
  124. 1300     'line between vertices m3 and m1 (case 3)
  125. 1310     X1=XH(M3): Y1=YH(M3): X2=XH(M1): Y2=YH(M1)
  126. 1320     GOTO 6000
  127. 1330 '
  128. 1400     'line between vertex m1 and side m2-m3 (case 4)
  129. 1410     X1=XH(M1): Y1=YH(M1)
  130. 1420     X2=(H(M3)*XH(M2)-H(M2)*XH(M3))/(H(M3)-H(M2))
  131. 1430     Y2=(H(M3)*YH(M2)-H(M2)*YH(M3))/(H(M3)-H(M2))
  132. 1440     GOTO 6000
  133. 1450 '
  134. 1500     'line between vertex m2 and side m3-m1 (case 5)   
  135. 1510     X1=XH(M2): Y1=YH(M2)
  136. 1520     X2=(H(M1)*XH(M3)-H(M3)*XH(M1))/(H(M1)-H(M3))
  137. 1530     Y2=(H(M1)*YH(M3)-H(M3)*YH(M1))/(H(M1)-H(M3))
  138. 1540     GOTO 6000
  139. 1550 '
  140. 1600     'line between vertex m3 and side m1-m2 (case 6)
  141. 1610     X1=XH(M3): Y1=YH(M3)
  142. 1620     X2=(H(M2)*XH(M1)-H(M1)*XH(M2))/(H(M2)-H(M1))
  143. 1630     Y2=(H(M2)*YH(M1)-H(M1)*YH(M2))/(H(M2)-H(M1))
  144. 1640     GOTO 6000
  145. 1650 '
  146. 1700     'line between sides m1-m2 and m2-m3 (case 7)
  147. 1710     X1=(H(M2)*XH(M1)-H(M1)*XH(M2))/(H(M2)-H(M1))
  148. 1720     Y1=(H(M2)*YH(M1)-H(M1)*YH(M2))/(H(M2)-H(M1))
  149. 1730     X2=(H(M3)*XH(M2)-H(M2)*XH(M3))/(H(M3)-H(M2))
  150. 1740     Y2=(H(M3)*YH(M2)-H(M2)*YH(M3))/(H(M3)-H(M2))
  151. 1750     GOTO 6000
  152. 1760 '
  153. 1800     'line between sides m2-m3 and m3-m1 (case 8)
  154. 1810     X1=(H(M3)*XH(M2)-H(M2)*XH(M3))/(H(M3)-H(M2))
  155. 1820     Y1=(H(M3)*YH(M2)-H(M2)*YH(M3))/(H(M3)-H(M2))
  156. 1830     X2=(H(M1)*XH(M3)-H(M3)*XH(M1))/(H(M1)-H(M3))
  157. 1840     Y2=(H(M1)*YH(M3)-H(M3)*YH(M1))/(H(M1)-H(M3))
  158. 1850     GOTO 6000
  159. 1860 '
  160. 1900     'line between sides m3-m1 and m1-m2 (case 9)
  161. 1910     X1=(H(M1)*XH(M3)-H(M3)*XH(M1))/(H(M1)-H(M3))
  162. 1920     Y1=(H(M1)*YH(M3)-H(M3)*YH(M1))/(H(M1)-H(M3))
  163. 1930     X2=(H(M2)*XH(M1)-H(M1)*XH(M2))/(H(M2)-H(M1))
  164. 1940     Y2=(H(M2)*YH(M1)-H(M1)*YH(M2))/(H(M2)-H(M1))
  165. 1950 '
  166. 2000 '
  167.  
  168. 6000 GOSUB 20000 'DrawIt -- vecout(x1,y1,x2,y2,z(k))
  169.  
  170. 7000 NEXT M 'case 0
  171.  
  172. 8000 'None in triangle
  173. 8001 NEXT K
  174.  
  175. 9000 'None in box
  176. 9001 NEXT I: NEXT J
  177. 9002 '
  178.  
  179. 10000 RETURN
  180.  
  181. 10999 '*****************************************************************
  182. 11000 'Subroutine
  183. 11010 H(M)=D(I+IM(M-1),J+JM(M-1))-Z(K)
  184. 11020 XH(M)=X(I+IM(M-1))
  185. 11030 YH(M)=Y(J+JM(M-1))
  186. 11040 RETURN
  187.  
  188. 12000 'Subroutine
  189. 12010 H(0)=(H(1)+H(2)+H(3)+H(4))/4   'average heights of four corners to get center height
  190. 12020 XH(0)=(X(I)+X(I+1))/2
  191. 12030 YH(0)=(Y(J)+Y(J+1))/2
  192. 12040 RETURN
  193.  
  194. 20000 'DrawIt subroutine 3D is squashed into 2D by ignoring z(k)-elevation
  195. 20010 LINE (X1,Y1)-(X2,Y2)
  196. 20020 RETURN
  197.  
  198. 30000 '*****************************************************************
  199.